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ABSTRACT 

Environmental two-level systems (TLS) have been identified as significant decoherence sources in Josephson junction (JJ) 
based circuits. For such quantum devices to be functional, the removal or control of the TLS is a necessity. Understanding the 
microscopic origins of the ‘strongly coupled’ TLS type is one current path of investigation to that end. The delocalized oxygen 
model suggests the atomic position of an oxygen atom is spatially delocalized in the oxide forming the JJ barrier. In this report 
we extend this model from its previous 2+1D construction to a complete 3D description using a Wick-rotated time-dependent 
Schrodinger equation to solve for time-independent solutions in three dimensions. We compute experimentally observable 
parameters for phase qubits and compare the results to the 2+1D framework. We devise a Voronoi classification scheme to 
investigate oxygen atoms delocalizing within strained and non-strained crystalline lattices, as well as realistic atomic positions a 
JJ amorphous tunnel barrier constructed in previous density functional studies. 


Introduction 

Decoherence sources such as the environmental two-level sys¬ 
tem (TLS)^’^ are presently a major hurdle to the realisation 
of operational superconducting qubits and other Josephson 
junction based quantum devices. Many properties of the TLS 
have been probed experimentally, demonstrating relatively 
long decoherence times whilst showing them to be stable and 
controllable.^"^ However, identifying their true microscopic 
nature remains an open problem. Controllable superconduct¬ 
ing qubit architectures (charge, flux and phase) make it possi¬ 
ble to study single ‘strongly coupled’ TLSs,^^’^’^ as opposed 
to weakly coupled ensembles of TLSs that may be responsible 
for 1// noise.Although the Standard Tunneling Model 
(STM)^^’ can be used to explain all existing experimental 
results, to improve devices we need a detailed microscopic 
understanding of these defects. 

Many microscopic models have now been proposed: de¬ 
localised oxygens,^^’^"^ surface state interactions,^^’po- 
laron dressed electronsand alien species (primarily hy¬ 
drogen) all of which are possible decoherence sources 
and/or explain the origin of strongly coupled TLSs. Better 
fabrication methods, higher vacuums and more robust qubit 
designs have historically suppressed or diminished the re¬ 
sponse of such noise sources.^^"^^ However, continuing down 
the path of a perfectly clean but amorphous dielectric may no 
longer be the optimal direction. Whilst crystalline layers are 
difficult to construct in this environment, they are possible 
and show a dramatic decrease (of up to 80%) in visible TLSs 
over their amorphous counterparts.^^ 

Such a decrease could be explained via our previous work 
on the delocalised oxygen model, where we suggest the 
origin of some TLS defects are oxygen atoms in a spatially 
delocalised state due to the non-crystalline nature of the dielec¬ 
tric layer. As such, this is a concrete microscopic ensemble 
of the classic atomic tunneling model: typically cited as moti¬ 


vation for the STM (which is more general, as the tunneling 
object need not be a single atom). 

Our previous model discussed a simplified configuration 
space using low dimensional arguments, culminating in a 
2+ID framework where a cage of six aluminium atoms sur¬ 
rounded a central oxygen atom in three dimensions and cal¬ 
culated observations based on the oxygen delocalising on a 
two dimensional plane. The dimensional simplicity of this 
model allowed a direct matrix diagonalisation approach to 
solve a time independent Schrodinger equation, which due 
to memory limits could not be extended into a complete 3D 
representation of a TLS defect. 

Here we introduce an alternative method for this model 
which can obtain a complete three dimensional picture of 
the delocalised oxygen ground and excited states. We use a 
Wick-rotated time-dependent Schrodinger equation to solve 
for time-independent solutions in three dimensions (see Meth¬ 
ods for more detail). As the work in this report is a higher 
dimensional extension of the 2+ID model, its terminology 
has been retained herein. It is suggested the reader be familiar 
with the discussions in Ref 14, although a short summary of 
the nomenclature is outlined below. Using this more powerful 
technique, we are able to move beyond a simplistic cubic 
representation and investigate true lattice configurations and 
realistic amorphous positions from our DFT studies. 

Two defect ‘types’ are investigated, both of which stem 
from deforming the crystalline lattice of corundum: the low 
temperature and pressure phase of AI 2 O 3 (the closest Al-O 
bond distance in this phase is 1.85 A). Defect type A; where 
the aluminium-oxygen bond distance is shortened, forcing 
the oxygen to occupy one of two off-axis positions, and defect 
type B; where the opposite occurs: the aluminium-oxygen 
bond distance is lengthened, allowing two preferred oxygen 
positions, on-axis. 

The trigonal lattice structure of this crystalline material 


and the amorphous nature of a JJ barrier give rise to many 
complications for a microscopic model, due to potential con¬ 
tributions from many neighbouring atoms to a possible defect 
site. Therefore, to fully understand the response to these in¬ 
teractions, we start with a very simplistic cluster of atoms 
and add complexity once an understanding of rudimentary 
behaviour is achieved. 

Defining the oxygen position to be at an origin, two in¬ 
troduced aluminium atoms can be considered as pairs (x — 
—X, -\-X = {±1.85 A}) lying on a cardinal axes; which are 
identified as |X| or simply referred to as the ‘defect pair’. Dis¬ 
placing \X\ equidistantly from this origin (i.e. moving away 
from optimal crystalline configuration) will yield either an A 
or B type defect, depending on the direction of displacement. 
\Y\ and |Z| describe aluminum pairs in the y and z directions 
respectively. 

Eigenenergies of a system are frequently discussed in a 
relative fashion using a convention where Etj = Ej — Ei, such 
that the ground state (Eq) to first excited state (Ei) energy 
splitting is defined as Eoi. 

Results 

The capability to investigate oxygen delocalisation in 3D, us¬ 
ing realistic atomic positions within an amorphous layer of 
AI 2 O 3 generates an extensive state space. Hence we start the 
investigation using a simple three atom system comprising 
of an oxygen atom and a confining aluminium pair: \X\. An 
eigenspectrum of the ground and five lowest excited energy 
levels of this system is generated by increasing or decreasing 
the pair separation distance from the corrundum lattice dis¬ 
tance 1.85 A which is depicted in Figure 1 (for comparison 
with the 2D model, see Figure 2 in Ref 14). 

In this figure an (an)harmonic division about |X| ~ 1.85 
A separates two regions with dissimilar properties. In the 
region where \X\ > 1.95 A, a degenerate Eqi is observed and 
as the separation distance has increased, can be identified as 
a possible B type defect. The second region (\X\ < 1.55 A) 
has a tightly bound eigenspectrum compared to the rest of the 
map although this yields no degeneracy. 

In three dimensions, the potential minima manifests as a 
sphere around the location of an atom. As a consequence the 
three atom Al-O-Al chain produces a unique, rotationally 
symmetric ground state which is reminiscent of an oxygen 
interstitial defect in crystalline germanium^^ and depicted 
in Figures 2a&b. Comparatively, this region in the 2D case 
indicated the presence of an A type defect, as the 2D potential 
minima ring would be projected onto a plane, collapsing a 
degree of freedom. 

Two other observations concerning the 2D to 3D transition 
is the dramatically reduced total energy values as wavefunc- 
tions are no longer artificially confined, and the higher excited 
states in the B type region now exhibit a quad degeneracy 
rather than two split doubly degenerate pairs - again because 
states are not artificially confined in the extra dimension. 



Figure 1 . Eigenspectrum of a three atom system Al-O-Al, 
over a varying distance separation \X\. Each excited state has 
been normalised with the ground state, which shows a clear 
Eqi degeneracy at large separation distance - indicative of a B 
type defect. An intermediate region exists as the separation 
distance decreases, which is approximately centered about 
the optimal Al -0 bond distance of corundum (1.85 A). At 
small separation distances, the eigenspectrum is distributed 
over comparatively small energy levels, but is not degenerate 
(see Figure 2). 

This simplistic three atom description can now be built 
upon to understand the interactions of a delocalised oxygen in 
a three dimensional space. A confining pair is introduced in 
an additional dimension to the configuration in Figures 2a&b, 
and the unique ground state is lifted to a degenerate Eqi pair; 
observing an A type TFS defect. Figures 2c&d illustrate this 
effect with a confining pair |F| = 3.5 A. This distance is an 
arbitrary choice, as the shift occurs even at the cut-off limits 
imposed by the potential portion of the model^^ (see Methods). 
In other words, any additional potential contribution which 
does not share the axial symmetry of the Al-O-Al chain will 
lift the degeneracy and result in a TFS. 

Oxygenic Orbitals 

A complete picture of the ground and excited state wave- 
functions in 3D is now possible, which allows us to further 
investigate the properties of the two-level system and illustrate 
the importance of crystalline dielectrics in future Josephson 
junction devices. 

A confined harmonic state can exist with many atomic con¬ 
figurations, and as with the low dimensional model, setting 
variables symmetrically is one of these cases. Figure 3 depicts 
the four lowest energy wavefunctions of an oxygen with six 
confining aluminium atoms: \X\ — \Y\ — \Z\ — 1.6 A. Here the 
oxygen atom has no additional local minima to occupy, iden¬ 
tifying this configuration as spatially localised, the harmonic 
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Figure 2. Ground (left) and first excited (right) state wavefunctions of an Al-O-Al chain with \X\ = 1.5 A (aluminium atoms 
are depicted as gray spheres). Top axes (a & b) show a unique ground state, the Eu pair is degenerate in this case. Two more 
aluminium atoms are introduced at \Y\ = 3.5 A on the bottom axes (c & d), which causes the degeneracy to lift to the Eqi pair. 


approximation holds and the system is not considered as a 
TLS candidate. 

From this localised case, we extend the \Y\ confinement out 
to 2.1 A. Using the Figure 1 eigenspectrum we can predict this 
configuration to be a B type defect; where the \Y\ separation 
distance is increased, and an oxygen dipole forms parallel to 
the y axis. This is indeed the case as shown in Figure 4. 

A particularly complex phenomenon emerged from the low 
dimensional model: quad degenerate ground states, where the 
energy difference between two degenerate pairs Eqi and E 23 
approach the difference of the pairs themselves, i.e. £'01 = 
£23 ~ Ei 2 . Figure 5 shows a configuration expressing this 
behaviour. 

The discussion in Ref 14 indicated that this form of degen¬ 
eracy has a low probability of occurrence compared to A or 
B type defects. Additionally, if Eu > 100 GHz the higher 
states (£2 and £ 3 ) can be ignored completely as energies of 
this magnitude dissociate Cooper-pairs, hence can be viewed 


as an operational upper bound for Josephson junction devices. 
Splitting energies in 3D are much smaller than the 2-i-lD case, 
which suggests that the probability of these defects being ex¬ 
perimentally visible may be higher than predicted with the low 
dimensional model based solely on this measure. However, 
a comprehensive analysis of the three dimensional configura¬ 
tion space would need to be undertaken before real estimates 
of this behaviour can be given. 

Another interpretation of this phenomena has to do with the 
energy resolution of the coherent single-pulse resonant driv¬ 
ing spectroscopy^ used to obtain qubit-TLS coupling (S^ax) 
values. The minimum resolvable energy using this process 
is approximately 1 MHz, which may indicate that Smax anti¬ 
crossing measurements are manifestations of qubit-TLS cou¬ 
plings at the £12 level, if a configuration such as Figure 5 is 
scrutinised. £01 and £23 splitting levels would be hidden in 
their degenerate ^bspace and manifest as an effective two- 
level defect with £01 =31.26 GHz. 
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Figure 3. Wavefunctions of an oxygen confined by six 
equidistant aluminium atoms at |X| = |F| = |Z| = 1.6 A. This 
configuration exhibits a unique ground state (as shown in the 
energy level diagram) and hence is considered to be spatially 
localised. 



|0) |1) |2) |3) 


E 23 = 36.37 kHz 
Ei 2 = 23.65 THz 
Eoi = 6.53 kHz 

Figure 4. Wavefunctions of an oxygen confined by six 
aluminium atoms, where the symmetry of Figure 3 has been 
broken. |X| = |Z| = 1.6, |F| = 2.1 A, which manifests as a B 
type defect. 

A type defects and dipoie considerations 

Whilst adding additional confinement pairs causes the unique, 
rotationally symmetric ground state of |X| < 1.55 A (Figures 
2 a&b) to become degenerate; the third dimension yields a 
complication in the dipole measurement for the A type region. 

Consider a system with parameters |X| = 1.53, |F| = 
1.97, |Z| = 1.95 A, illustrated in Figure 6. This system ex¬ 
hibits TLS behaviour, with £'01 =40.63 MHz, and a dipole 
strength in y of 0.26 eA: perpendicular to the confining x axis 
as expected. The leftmost axis of Figure 6 shows a 2D projec¬ 
tion of the first excited state, illustrating no major differences 
in the response of the 3D and 2+ID models. 

However, a small change in the y confinement alters the 
system in a non-trivial manner. Moving |F| from 1.97 to 1.9 
A (for example) crosses a bifurcation in state space. As |Z| 
is now the least confining pair, the dominant dipole direction 
flips to z as shown in Figure 7. 

Ultimately this does not add any complexity to the model: 


Figure 5. Wavefunctions of an oxygen confined by six 
aluminium atoms, where the confinement configuration is 
studiously pathological. \X\ = 1.6, |F| = 2.1, |Z| = 2.0 A, 
yielding a quad degeneracy in the ground state due to the £12 
splitting existing below the superconducting gap (~ 100 GHz, 
see text). 

minimal confinement in y generates an A type defect with 
a dipole in y, perpendicular to x; minimal confinement in z 
generates an A type defect with a dipole in z, perpendicular to 
X. As strongly coupled defects are experimentally identified 
via avoided level crossings in qubit spectroscopic diagrams,^ 
the dipole moment of the defect couples to the electric field 
across the junction.Therefore, to be identified as a TLS, 
the defect must be aligned to this field. It follows then, that 
the dipole directions for an A type defect in this model can be 
considered equivalent. Additionally, as with B type and any 
other possible TLS defect, both A type alignments will not 
appear on a spectroscopic scan as an avoided level crossing if 
their dipole is perpendicular to the external field. 

Defects in corundum 

Results such as Figures 1 and 3 identify configurations which 
exhibit harmonic eigenspectrums and localised oxygenic po¬ 
sitions. This behaviour is expected when an oxygen atom is 
confined appropriately and has no reason to tunnel. One way 
to check this assumption (and consequently further validate 
the TLS model) is to calculate the potential landscape seen by 
a single oxygen atom from bulk corundum (a-Al 203 ). The 
TLS model, delocalising an oxygen across all potential space 
about its minimum energy position, should yield a localised, 
harmonic wavefunction positioned at the location where the 
oxygen was removed (Le. the lattice position associated with 
the vacancy). 

An oxygen atom in a corundum lattice is identified by a 
Voronoi polytope (see Methods) and a nearest neighbour atom 
cluster is used to calculate the wavefunctions and energies 
of an oxygen atom located inside the potential landscape 
generated by the corundum lattice. As predicted, the result is 
a localised oxygen atom locked into the lattice position of the 
vacancy, and the oxygenic orbitals behave in a similar fashion 
to the harmonic system calculated previously in Figure 3 (data 
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Figure 6. Wavefunctions of an A type defect with 
confinement \X\ = 1.53, \Y\ = 1.97, |Z| = 1.95 A which 
yields dominant dipole in the y direction with a strength of 
0.26 eA. Leftmost axis shows an xy projection of the first 
excited state to compare with 2+ID results in Ref 14. 

not shown). 

Straining the Locai Lattice 

In an attempt to put bounds on the amount of amorphousness 
required for a near-crystalline structure to exhibit TLS prop¬ 
erties, we apply a localised strain tensor on the defect cluster 
and observe its effects. Two strain ratios ^ = ^0.019 A rep¬ 
resenting a compressive force along the x axis, and a tensile 
force also along x respectively, attempt to simulate a local 
strain from the localised Al-0 separation distance of 1.85 A 
in the direction of the A type and B type regions respectively 
(see Figure 1). 

The original crystalline system has a ground to first excited 
state splitting of £’oi = 13.48 THz. Comparing the £”01 values 
of the strained results (compressive: 35.42 THz, tensile: 34.84 
THz) to the crystal result, we see an increase in splitting 
energy (data not shown). This response is opposite to the 
simple straining of Figure 1 where the £01 energy was a 
maximum at the crystalline Al-O separation distance. This 
suggests that local Poissonian strain does not give rise to 
immediate TLS behaviour—at least along the sampled axis x. 

One possible explanation of this result comes from our pre¬ 
vious study concerning strain on the 2+ID model: there 
exists a preferred strain axis which responds rapidly in com¬ 
parison to other directions, although this axis is difficult to 
acquire directly from a completely harmonic case. It is there¬ 
fore pertinent to investigate the amorphous configurations in 
more detail. 

Defects in Josephson Junctions 

Applying the Voronoi classification scheme the myriad atomic 
configurations within the Josephson junction models con¬ 
structed in Ref. 27 allows us to investigate realistic amorphous 
clusters residing within a JJ tunnel barrier. 



Figure 7. Wavefunctions of an A type defect with 
confinement \X\ = 1.53, \Y\ = 1.9, |Z| = 1.95 A which yields 
dominant dipole in the z direction with a strength of 0.21 eA. 
Leftmost axis shows an xz projection of the first excited state 
to compare with 2+ID results in Ref 14. 

Using the structure with an AIO 1.25 amorphous barrier with 
a density 0.75 that of corundum we classify clusters for all 
oxygen atoms extant in the tunnel barrier. It is posited that a 
large void space around an oxygen may allow it to delocalise 
over a large region and manifest as a B type defect (however, 
this is pure speculation). The volume of each Voronoi cell 
can be calculated (as the cell is essentially a convex hull of 
some vertices describing a Delaunay tessellation) and we 
choose the five cells with the largest volume to investigate. 
Another conjecture is that an asymmetric cell (obtained by 
some form of sheer in the local cluster configuration) yields 
TLS behaviour, thus a final cluster is investigated which has 
the most asymmetry - identified using the length parameter 
of the spherocylinder formalism outlined in 30. The five cells 
with the most asynunetry were also identified; these clusters 
and their properties are displayed in Figure 8 along with the 
large volume candidates. 

As with the cluster results of the crystalline and Poisson 
strain investigations, the few amorphous clusters studied failed 
to identify any TLS behaviour. Whilst each oxygen displays 
a different spectra, none of the £01 splitting are degenerate 
and dipole strengths are no larger than 0.034 eA (from the red 
cell, identified as the oxygen with the most asymmetry in its 
cluster) in any cardinal direction. 

Discussion 

Using the direct diagonalisation method in three dimensions 
is seennngly intractable with current memory limitations, al¬ 
though this does not exclude investigations in the 3D do¬ 
main. The Wick-rotated time-dependent Schrodinger equation 
method has proved to be a valid tool to study this space, albeit 
with computational time disadvantages. 

Whilst the 2+ID model allows us to investigate phase space 
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Figure 8. Top: Oxygen atoms with the 5 largest Voronoi 
volumes and the 5 most asymmetric Voronoi cells are 
highlighted in within a Josephson junction comprised of 
aluminium (gray) and oxygen (red) atoms. Bottom: 
Eigenspectrum of the ground state and four lowest excited 
state values for each of the candidate oxygen atoms. 


The free parameter concern stated in Ref. 14 is somewhat 
mitigated through this process, as clusters with many atoms 
can be described with one Voronoi cell (clusters of up to 24 
atoms are described in Figure 8 for example). 

The affect of local strain on a crystalline cluster observed 
behaviour which does not align with the simple strain inves¬ 
tigations undertaken with the idealised 6 atom model, sug¬ 
gesting that the deformation axis chosen is not the preferred 
strain axis identified in Ref. 14 and that localised strain may 
not lead to TLS behaviour directly. 

With such a small sample size it is difficult to draw any 
direct conclusions from the results in Figure 8. Other classifi¬ 
cation metrics or methods of identification such as extending 
the Euclidian distance Voronoi classification discussed here 
to the Voronoi S-networks considered in 30 may yield TFS 
behaviour with the present model. 

However, the more probable explanation is that the sur¬ 
rounding atom clusters exert enough confinement energy to 
the oxygen atoms to keep them spatially localised, unlike the 
simpler 6 atom confinements. One phenomena which may as¬ 
sist in successful TLS classification, is the inclusion of many 
body interactions (the C type defect in Figure 1 of Ref. 14). 
As oxygen atoms tunnel through potential barriers, the sur¬ 
rounding atomic neighbours may deform and create a multi¬ 
atom TLS candidate. Such behaviour has been studied pre¬ 
viously,^ and multi-body identification approaches have 
recently been investigated.^^ Generalising the delocalised 
oxygen atom to a multi-body description is a topic for future 
work. 


Methods 

To identify A and B defect types in either amorphous or 
strained crystalline systems, an effective single particle Hamil¬ 
tonian can be used. This simplifies the system by ignoring 
any time evolution properties or many-body interactions and 
yields the equation 


in a much more flexible way, one of the major variations in the 
models is the final energy values, as the 3D model is no longer 
collapsed into fewer dimensions. This result suggests that if 
one is attempting to identify exact values for a delocalised 
oxygen TLS, they should be utilizing the 3D method. Con- 
trarily, if trends are required, the 2-t-lD model provides the 
more efficient solution. Whilst splitting energies and dipole 
results may not completely reflect the atomic positions of the 
surrounding clusters in this case, it is evident from discussions 
like those concerning Figures 6 and 7 that there are 2D to 3D 
equivalencies. 

Using a Voronoi classification scheme to identify atom 
clusters which describes the local potential environment of 
an oxygen atom is shown to be a powerful method, describ¬ 
ing a lattice position in crystalline corundum as well as the 
properties of any selected oxygen atom in a model junction. 



^^oxy 


+ nr), 


( 1 ) 


where moxy is the mass of an oxygen atom and V(r) is 
the potential due to the surrounding lattice. This potential 
is modelled by the empirical Streitz-Mintmire potential,^^’ 
which adequately captures the variable oxygen charge states 
occurring when oxygen is present in a predominantly metallic 
environment like a JJ metal-oxide interface. A complete de¬ 
scription of the parameters used can be found in Ref 14 as the 
implementation remains consistent with this model. 

To obtain full 3D eigenstates of this Hamiltonian we imple¬ 
ment a Wick-rotated time-dependent Schrodinger equation to 
imaginary time t = 
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(3) 


which yields a general solution to the wavefunctions 


'i>{r,T)=Y,anWn{r)e (4) 

n=0 

Here, Un are coefficients based on initial conditions of the 
system where n — 0 indicates the ground state, n = 1 the first 
excited state etc. and En is the corresponding eigenenergy. 
As £”0 < £”1 < £2 < ..evolving Eq. 4 to large imaginary 
time will provide a good approximation to the ground state 
influenced by the time-independent potential V (r). 

This solution is obtained by numerically approximating 
the spatial derivatives with finite differences. For stability 
and precision, we descritise the space using a 7-point central 
difference method 


order screening of the surrounding lattice potential, represent¬ 
ing an acceptable approximation to higher order screenings 
(e.g. also including second nearest neighbours, which add 
appreciable computational intensity). Clusters identified by 
this method use the Streitz Mintmire potentials generated by 
GULP,"^^ which also account for monopole-monopole and 
self energy interactions to counteract any net charge generated 
by selecting an arbitrary cluster of atoms from a lattice with 
periodic boundaries. 

Arbitrarily straining a corundum crystal, which has trigonal 
symmetry and six elastic constants,'^^ will exert force on the 
structure in a non-trivial manner and act on many different 
crystal planes. The materials science complications of this 
phenomena are well beyond the scope of this investigation, 
and instead we approximate a local strain by the Poisson ef- 
fect"^^ using the values of the adiabatic bulk and shear moduli 
to obtain a Poisson’s ratio v = 0.234. This value is calcu¬ 
lated as an equivalent isotropic polycrystalline aggregate to 
Sapphire (a mineral variety of corundum). 

The relationship 



(7) 


/"(ro) » 

2/_3 - 27/_2 + 270/_i - 49O/0 + 270/1 - 27/2 + 2/3 


where fk — f (ro + kh), calculated with a step size /i = 0.01 

A. 

Excited states are obtained by Gram-Schmidt orthogonali- 
sation:^^’'^^ choosing eigenfunctions within a degenerate sub¬ 
space to be orthogonal to the system’s basis. For example, 
consider a converged, orthonormal ground state |i//b) and an 
initial, non-orthonormal guess for the first excited state | ). 

Subtracting the projection of the excited state along the ground 
state from the initial guess yields a vector orthogonal to the 
ground state 

\^i} = \v[)-\Vo}{Vo\v'i)- ( 6 ) 

This resultant vector can then be normalised and evolved 
through the same process as the ground state until convergence 
is achieved. Similarly, the second excited state requires the 
same treatment, although projections along both |i//b) and 
|i/ri) are needed. 

An active repository of the software which implements 
these methods can be found via Ref 41. 

When dealing with the corundum crystal and amorphous JJ 
barrier clusters, we apply a Voronoi tessellation^^ to the lattice 
coordinates of each system using a Euclidean distance metric. 
A polytope representing the convex hull encompassing an ori¬ 
gin position symbolises the oxygen vacancy. Atoms situated 
within polytopes sharing an edge with the origin polytope are 
selected as the model cluster. These atoms constitute a first 


can be used to strain our atom cluster for example, with a 
length increase of AL in the x direction, and a length decrease 
of AL' in the y and z directions. 
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